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The static and dynamic structure of liquid Al is studied using the orbital free ab-initio molecular 
dynamics method. Two thermodynamic states along the coexistence line are considered, namely 
T = 943 K and 1323 K for which X-ray and neutron scattering data are available. A new kinetic 
energy functional, which fulfills a number of physically relevant conditions is employed, along with a 
local first principles pseudopotential. In addition to a comparison with experiment, we also compare 
our ab-initio results with those obtained from conventional molecular dynamics simulations using 
effective interionic pair potentials derived from second order pseudopotential perturbation theory. 



I. INTRODUCTION. 

Molecular dynamics (MD) methods have a long tradi- 
tion as a useful technique to study the properties of liquid 
systems, and the last fifteen years have witnessed a large 
spread on the application of ab-initio molecular dynamics 
methods, based on the density functional theory (DFT). 
This theory allows calculation of the ground state elec- 
tronic energy of a collection of atoms, for given nuclear 
positions jd'El and also yields the forces on the nuclei via 
the Hellmann-Feynman theorem. It enables to perform 
molecular dynamics simulations where the nuclear posi- 
tions evolve according to classical mechanics whereas the 
electronic subsystem follows adiabatically. 

In this paper we present the results of an ab-initio 
molecular dynamics simulation on the static and dy- 
namic properties of liquid Al at thermodynamic condi- 
tions around the triple point. Liquid aluminium has usu- 
ally been considered as a simple metal in which the core 
electrons forming the ion can be clearly distinguished 
from the valence electrons, and moreover the core elec- 
trons do not significantly overlap with those of neigh- 
bouring ions. Therefore, the system consists of a binary 
mixture of ions and valence electrons, where the former 
may be treated classically whereas the electrons must be 
treated quantum mechanically. 

However, in wide regions of the density-temperature 
plane, the simple metals have usually been treated as 
an effective one-component fiuid of ions interacting by 
means of density-dependent effective interionic pair po- 
tentials, derived from ionic pseudopotentials by applying 
second order perturbation theory. This approach, which 
will be refered to as the linear response theory (LRT) , has 
often been used as the starting point for the study i 
static and dynamic properties of the simple metals.l 
It has also been the approach followed in most studies on 



* e-mail: david@liql.fam.cie.uva.es 



the static structure of liquid aluminium. BBilli Among 
them, we mention the work of Dagens et aln who ob- 
tained an effective interionic pair potential, derived from 
a non-local pseudopotential which was constructed from 
the valence charge density induced by an Al"'"'^ ion placed 
in an electron gas at the metallic density. From this 
potential, Jacucci et aln calculated the static structure 
factor of liquid aluminium by means of MD simulations; 
their results showed fair agreement with experiment, with, 
a main peak somewhat higher. Also Hafner and Janktll 
have studied the liquid static structure of aluminium by 
means of an effective interionic pair potential derived 
from an (A^initio pseudopotential originally developed by 
Harrison|13 whereas the corresponding liquid static struc- 
ture was derived by means of MD simulations. In their 
calculation of the pseudopotential, the authors used the 
coefficient for the exchange-correlation potential between 
the core and valence electrons as a fitting parameter in 
order to obtain agreement with the experimental static 
structure factor. 

Whereas the previous work dealt only with the static 
properties of liquid Al near melting, the work by Ebbsjo 
et aln also considered some dynamic properties. In fact, 
these authors performed MD simulations for three differ- 
ent interionic pair potentials, two of them based on non- 
local pseudopotentials and the other one based on the 
local Ashcroft's pseudopotential,ll3 which showed rather 
different shape, specially outside the repulsive core. De- 
spite those differences, all them gave fairly similar results 
for the liquid static structure which agreed well with the 
experimental data whereas the main discrepancies ap- 
peared in the dynamic structure. 

A rather different aprirpach has been followed by Chi- 
hara and coworkers OtJ Their quantum hypernetted 
chain (QHNC) method treats the ions and electrons on 
basically equal footing by combining liquid state integral 
equations with the density functional formalism. More- 
over, it does not rely on the pseudopotential ideas, gives 
rise to a self-consistent scheme to determine the liquid 
static structure and yields an effective interionic pair 
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potential which depends on the particular liquid static 
structure. 

Although the LRT approach has produced reasonable 
results for the liquid alkali metals, when the valence of 
the system is increased, its validity becomes more ques- 
tionable. In addition, even for the alkali metals the LRT 
is less justifiable for thermodynamic states approaching 
the critical point, and it is certainly wrong near the crit- 
ical point. This limitation of the "standard" theory has 
estimulated the j.isn of fkat-principles molecular dynam- 
ics techniques ,E3EZlliJO'E£l where the electronic density, 
total energy and forces are obtained by using the Kohn- 
Sham (KS) formulation of the density functional theory 
(DFT).cl However, the computational demands of these 
ab-initio methods, where KS orbitals are used to describe 
the electronic density and to compute exactly the elec- 
tronic kinetic energy, grow very rapidly with system size, 
and their memory requirement is also quite large. These 
considerations have restricted the sizes of the systems 
studied so far, to about 60 atoms, and have limited sim- 
ulation times tn_axQiiad_2-5 picoseconds in the cases of 
Rb, Cs and Hg,E30ll3'El and 64-128 atoms, with simu, 
lation times of 0.15-0.85 picoseconds in the case of Na.E£l 
These limitations can be at least partly overcome if the 
exact calculation of the electronic kinetic energy is given 
up in favour of an approximate kinetic energji-functional 
of the electronic density. Within this schemeEiJ the num- 
ber of variables describing the electronic state are enor- 
mously reduced, especially for large systems, enabling the 
study of larger systems for longer simulation times. This 
approach has already been used for several studies on 
solids,c3 clusters pAfiTid some liquid metals (Li, Na_^g, 
Al) near melting.oE^ We have recently presentedcZl an 
application of this method to study the static structure 
and some dynamic properties of expaaded liquid Cs, for 
which experimental data is available.E3 In that study 125 
particles were used and the simulation time was 17-35 
picoseconds after an equilibration time of 11-25 picosec- 
onds. Also, another studycS with 205 particles for liquid 
Al gave results for the static structure in|-gpod agreement 
with experiment. Recently, Anta et alEB have also ap- 
plied the same scheme to study the ionic and electronic 
static structure of liquid Al near melting, leading to re- 
sults for the static structure factor in excellent agreement 
with experiment. 

The static structure factor of liquid A1 .has been mea- 
sured by both neutronoEJ and X-raytirEj diffraction. 
The dynamical structure of liquid Al near the triple point 
has also been investigated recently by Scopigno et al.L3 
using inelastic X-ray scattering (IXS). Note that the high 
value of the adiabatic sound speed for liquid Al ( « 4800 
m/s), prevents the use of the inelastic neutron scatter- 
ing (INS) technique for investigating the collective ex- 
citations for small values (roughly, for q < Qp, with 
qp « 2.70 being the main peak position of the static 
structure factor). Those IXS experiments have investi- 
gated the wavevector region 0.05gp < q < 0.5gp, obtain- 
ing several dynamical features previously observed in the 



liquid alkali metals, such as the existence of collective ex- 
citations up to (j-values larger that 0.5qp, which exhibit 
a positive dispersion in the sound velocity with respect 
to the hydrodynamic value. 

The layout of the paper is as follows. In section II 
we briefly describe the theory used in the orbital-free 
ab-initio molecular dynamics (OF-AIMD) simulations, 
giving some technical details, and focusing on the two 
problematic issues, namely, the kinetic energy functional 
and the local pseudopotentials needed to characterize the 
ion-electron interaction. In section III we present and 



discuss the results of the ab-initio simulations; moreover 
they are compared with further classical molecular dy- 
namics (CMD) simulations that we have performed based 
on LRT and the QHNC potentials, and with the avail- 
able experimental data. Finally some conclusions are 
drawn and possible ideas for further improvements are 
suggested. 



II. THEORY. 

The total potential energy of a system of N classi- 
cal ions enclosed in a volume V, and interacting with 
iVe — NZ valence electrons through a local electron-ion 
potential v{r), is written, within the Born-Oppenheimer 
approximation, as the sum of the direct ion-ion coulom- 
bic interaction energy, and the ground state energy of 
the electronic system subject to the external potential 
created by the ions, Vcxt{r, {Ri}) = YliLi v{\r- Ri\) , 



E{{Ri}) = J2 



\Ri — Ri 



+ EMr),V,^,{r,{Ri})], 



(1) 



where pg {r) is the ground state electronic density and Ri 
are the ionic positions. According to LRT, the ground 
state electronic density is given, in reciprocal space, by 



„LRT. ^ ^ 



n-^-{q)^F{q)n^^^{q) (2) 



n^^^^(g)-x(9,Po)i'(9) 



(3) 



where ^(g, po) is the response function of a uniform elec- 
tron gas of density po = NZ/V . Accordingly, the ground 
state electronic density is a superposition of spherically 
symmetric pseudoatomic densities around each ion, i.e.. 
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,LRT 



r-R, 



(4) 



and the electronic ground state energy is 
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0ind(g) = x{q,Po)v^{q) 



(5) 
(6) 



where Ey[po] is a structure-independent term. Within 
the LRT, the total potential energy can be written as 
a sum of a structure-independent term and a sum over 
pairs of an effective interionic pair potential 0off(^) = 

^V-R+0ind(i?)- 

Alternatively, DFT shows that the ground state elec- 
tronic density can be obtained by minimizing the energy 
functional E[p\, and the minimum value of the functional 
gives the ground state energy of the electronic system. 
The energy functional can be written 



E[p{r)] = [p] + E,^t [p] + Eh [p] + E^, [p] (7) 

where the terms represent, respectively, the electronic ki- 
netic energy, [p] , of a non- interacting system of density 
p^r) , the energy of interaction with the external potential 
due to the ions, 



Eo^t[p] = J drp(f)14xt(r), (8) 
the classical electrostatic energy (Hartrec term), 



Eh[p] = ^J Jdrds^^, 



(9) 



and the exchange-correlation energy, i?xc [p] , for which we 
will adopt the local density approximation. 



A. Technical details 

Given an explicit functional Ts[p], we can proceed to 
minimize Eg with respect to the p(r), but in order to 
maintain p{r) > everywhere, we have used as our 
system variable, an effective orbital, ip{r), defined as 
p(r) = ipir)'^^ with real We expand V'(^) in plane 

waves compatible with the simple cubic periodic bound- 
ary conditions of the simulation: 



G 



1 



G = ^ (711,712,713) . 



(10) 

(11) 
(12) 



where L stands for the side of the cube. This expansion 
is truncated at wavevectors corresponding to a given cut- 
off energy, i?cut , whose value is given in Table |. A real 



■0 implies that c_q = c^, with a real cq; consequently 
only the half-set {c^}'s need be treated as variables. 

The energy functional must be minimized with the nor- 
malization constraint t/[p(r)] = ^ydrp(f) = N^. which is 
imposed via the Lagrange multiplier /i, leading to the 
Euler-Lagrange equation 



_ S[E - pG] 



5p{r) 



SE 



5p(f) ^pi'P) 



p = p{r)-p = Q (13) 



for the ground state density. The minimization is per- 
formed with respect to the {c^jj's, instead of the elec- 
tronic density, leading to the equations: 



f drp{f)^{r) - 2pVc^ = 

OCo Jv 



dT 



4 / dfp{r)iP{r)e''^-^ - ApVc^ = (14) 



for the ground state density. The minimization of the 
functional is performed every time step of the simula- 
tion, using a simple quenching method: a fictitious "co- 
efficients' kinetic energy", T = ^Afc X^g I'^gI^' intro- 
duced, where Mc is the "coefficients' mass", and the dot 
denotes the derivative with respect to the fictitious "co- 
efficients' time", tc- This kinetic energy, rewritten in 
terms of the set {c,^}, together with the "potential en- 
ergy" JF, leads to the following "equations of motion" 
(Vcg e {eg}) 



-2 / df p[r)il){r)e 
h 



iG-r 



2pVc^ 



(15) 



These equations are-solved numerically using the Ver- 
let leapfrog algorithms with an electronic timestep A^c- 
The velocities are quenched at every step until the min- 
imum is reached within preset tolerances on T and the 
gradient of J-. The chemical potential p is not known in 
advance of the minimization, but replacing p in eqn. ( p^ ) 
by its stationary value J d'rp{r)n{r) / J dfn(r) at each 
timestep, gives good convergence to the ground state. 
For the present simulations, we have used Mc = 1.85 x 10^ 
hartree x (a.u.)'^ and a Ate = 1 x 10~* ps. 

The interatomic forces are obtained from the elec- 
tronic ground state via the Hellman-Feynman theorem, 
F, = -~Vj^^Eg[p{r),{Ri}], { i = 1---N) and Newton's 

equations, d^Ri/dt^ = Fi/Mi, are solved numerically for 
the motion of the ions using the Verlet leapfrog algorithm 
with a timestep At = 1.5 xlO^^ ps. 



B. The kinetic energy functional 

The kinetic energy functional, T^, is a critical ingredi-, 
ent of the energy functional. It is generally consideredE^ 
that the von Weizsacker term, 
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Twm] = ^ / dr\Vp{r)\ypir}, 



(16) 



is essential for a good description of the kinetic energy. It 
applies in the case of rapidly varying densities, and it is 
exact for one or two-electron systems. Further terms are 
usually added to the functional in order to reproduce cor- 
rectly some exactly known limits. In the uniform density 
limit, the exact kinetic energy is given by the Thomas- 
Fermi functional. 



TTF[pir)] 



3 

10 



dr p{r)kp{r) , 



(17) 



where kpir) — (37r^)^/'^(0(r)^/'^ is the local Fermi 
wavevector. In the limit of almost uniform density, LRT 
is correct, with a response function corresponding to a 
non-interacting uniform electron gas, given by the Lind- 
hard function, Xl{<1,Po)- 

Stimulated by the advantages of the orbital-free ab- 
initio simulations, there has been a renewed interest 
in the development of accurate kinetic eoergy function- 
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als. With_E 

coworkersEZtS have developped functionals which cor- 
rectly-j-ecover the Thomas-Fermi and linear response 
limita£!l and have also includedjthe quadratic response. E3 
Later, Carter and coworkera^J investigated these func- 
tionals and proposed a linear combination of them as a 
suitable form for T^; more recently they have also de- 
rived another expression which includes density depen- 
dent kernels. Unfortunately, an undesirable feature of 
these functionals is that they are not positive definite, so 
that minimization of the energy functional can lead to an 
unphysical negative kinetic energy. _ 

Chacon, Alvarellos and TarazonacJ have developed a 
different type of kinetic energy functional, which employs 
an "averaged density" and recovers the uniform and LRT 
limits. Their functional has beep_investigated and gener- 
alized by Garcfa-Gonzalez et aZ.EJ These functionals have 
the merit of being positive definite, but they are some- 
what complicated to apply and require order N more 
Fast Fourier Transforms (FFT's) than simpler function- 
als, and this diminishes the advantage of the orbital free 
approach over the full Kohn-Sham method. 

In this paper we. use a simplification of the averaged 
density approach,cil with the kinetic energy given by 





T,=Tw[p]+Tp[p] 


(18) 




Tp[p] dfp{r)'/^-'^k{r)^ 


(19) 


fc(f) = 


= {2k°pf j dsk{s)wp{2k^p\r~s\) 


(20) 




k{r) = (iTT^f'^p{r)P. 


(21) 


where kp is 


the Fermi wavevector corresponding 


to a 



mean electron density pq, and wp{x) is a weighting func- 
tion, determined by requiring the correct recovery of the 



LRT and uniform density limits. Note that k{r) appears 
as a convolution which can be performed rapidly by the 
usual FFT techniques. This functional is a generaliza- 
tion of one with (3 =j-J/3, used earlier by us in a study of 
expanded liquid CsBJ 

The details of the functional are given in appendix 
and its main characteristics are as follows: (i) /3 is a 
real positive number whose maximum value leading to 
a mathematically well behaved weight function is « 0.6, 
(ii) the functional recovers the uniform and LRT limits, 
and is positive definite, (iii) when fc^ — > because the 
mean electron density vanishes, e.g. for a finite system, 
the von Weizsacker term is recovered if (3 = 4/9, whereas 
for other values of /3, the limit is Tw -f CTtF: (iv) for 
values of /3 > 0.5 it is expected that p{r)ip{'r), which 
is the driving force for the dynamic minimization of the 
energy, (see eqn. ( p^ ) remains finite even for very small 
electronic densities p{r). 

Those two last properties will be important in the 
case of expanded liquid metals because of the appare- 
ance of large inhomogeneities in the atomic distribution, 
and therefore in the electron density, with regions where 
it becomes very small. Indeed, this situation has already 
been obsewed in the ab-initio simulations of expanded 
liquid NaB3 In systems for which the appearance of iso- 
lated atoms or clusters is likely the von Weizsacker term 
would be appropiate, and a functional with a value of /3 
as close as possible to 4/9 would be recommended. 

In the present simulations we have used P = 0.51, 
which in the limit pQ gives C= 0.046 and guarantees, 
at least for the thermodynamic states considered, that 
fi{r)ip{'r) remains finite and not too large everywhere so 
that the energy minimization can be achieved. 



C. Pseudopotentials 

Ab initio simulations using the full Kohn-Sham ap- 
proach (KS-AIMD) usually employ nonlocal pseudopo- 
tentiaJa. obtained by fitting to properties of the free 
atom.cj In an orbital-free approach where the electronic 
density is the variable, such non-local pseudopotentials, 
which act differently on different angular momentum 
components of the orbitals, cannot be used. Instead, 
local pseudopotentials must be developed which include 
an accurate description of the electronic structure in the 
physical circumstances of interest. 

When constructing a pseudopotential to be used for a 
liquid metal, it seems more appropriate to use a reference 
state which closely resembles the environment of an atom 
in the metal, which is quite different from free space. 
The pseudopotential used in this simulation has been 
obtained using the neutral pseudoatom (NPA) mcthodcj 
in which the reference state is an atom at the centre of 
a spherical cavity in the positive background of a uni- 
form electron gas. The density of the gas is taken to 
be the mean valence electron density of the system of 
interest, in our case, the liquid metal in a specific ther- 
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modynamic state. The radius of the cavity is such that 
the total positive charge removed from the hole is equal 
to the valence of the atom. First, a full Kohn-Sham den- 
sity functional (KS-DFT) calculation is performed to ob- 
tain the displaced valence electron density, rtps(r), i.e. 
the change in the electron density induced by the atom 
and the cavity. After pseudizing the nps{r) by eliminat- 
ing the core-orthogonality oscillations, an effective local 
pseudopotential is constructed which, when inserted into 
the uniform electron gas along with the cavity, repro- 
duces the displaced valence electron density previously 
obtained. 

Two approaches will be followed in this paper. The 
first uses LRT to reproduce the displaced valence elec- 
tron density leading to a LRT-based local pseudopoten- 
tial (LRT-PS) from which an effective interatomic pair 
potential is derived (see eqn.(||)) to be used in CMD sim- 
ulations; for further details we refer to Ref. The 
second approach uses the orbital-free density funtional 
theory (OF-DFT) to reproduce the displaced electron 
density and it is suited for OF-AIMD simulations. The 
development proceeds as follows. When the functional 
derivatives of the energy functional are performed, the 
Euler equation, eqn. (|l3|), for our pseudopotential in the 
jellium- vacancy system becomes 

^^sir) + V^^tir) + Vnir) + V^r) - ^l = , (22) 

where each of the terms is the derivative of the corre- 
sponding term in eqn. (|^), namely. 



vapor coexistence line ( 943 K and 1323 K), for |55diich_X=. 
ray and neutron diffraction data are available .cJCjly'EHl 
Table | gives further details on the thermodynamic states 
and several simulation parameters. In addition, we have 
also carried out classical MD simulations, using effec- 
tive interionic pair potentials derived from standard sec- 
ond order pseudopotential perturbation theory, with the 
LRT-PS's constructed as previously described (see also 
Ref. ^3|) and with pair potentials derived from the 
QHNC. 

In the OF-AIMD simulations 500 particles were 
treated in a cubic cell of the size appropriate to the den- 
sity, whereas more particles were used for the CMD sim- 
ulations (see Table |). In both sets of simulations, liquid 
static properties were evaluated (pair distribution func- 
tions and static structure factors) and several dynamic 
properties, both single-particle ones (velocity autocorre- 
lation function, mean square displacement) and collective 
ones (intermediate scattering functions, dynamic struc- 
ture factors, longitudinal and transverse currents). The 
calculation of the collective dynamic properties required 
long simulation runs in order to accumulate reasonable 
statistics; for example the OF-AIMD simulations lasted 
for 2 •10'' steps which correspond to 30 ps of simulation 
time. On the other hand, the CMD simulations run for 
10^ steps, amounting to 200 ps. 



(23) 



A. Pseudopotentials. 



with the expressions for the von-Weizsacker term and the 
/3-term given in appendix 

Vc^tir) = Vps{r) + Vcav{r) + Vjc\i{r) , (24) 



Vnir) = / ds>(s)/|f-s1, (25) 



with p(r) = po + n{r), and V^cij") is the exchange- 
correlation potential, obtained from the functional 
derivative of E^c[p\ evaluated at p{r). 

Due to the spherical symmetry of the system all the 
magnitudes depend just on r. Given p{r), Wps(?') can be 
obtained from eqn. (E4), and the constant /i is just an en- 
ergy origin which is set so as to obtain a pseudopotential 
that decays to zero for large distances. The pseudopo- 
tential constructed in this way, will be refered to as the 
OF-DFT-based pseudopotential (OFDFT-PS). 



III. RESULTS AND DISCUSSION. 

We have performed OF-AIMD simulations for liquid Al 
at two different thermodynamic states along the liquid- 



The local pseudopotentials described in section II C 



were constructed using a reference system mimicking the 
complex system to be studied. The pseudopotentials 
change with the thermodynamic state considered and 
therefore are not transferable to other states. Figure |l| 
shows the Fourier transforms of the non-coulombic part 
of the pseudopotentials obtained from the LRT and OF- 
DFT approaches outlined above. The two schemes lead 
to similar pseudopotentials with the main differences be- 
ing at low g-values and in the amplitude of the oscilla- 
tions at large q. Note that in both approaches the same 
pseudized displaced valence electronic density of an atom 
in a jellium- vacancy model is reproduced, although the 
OF-DFT was used in one case and the LRT in the other. 
Consequently the differences in the two pseudopotentials 
should reflect the importance of nonlinear effects which, 
according to the present results seem to be more impor- 
tant at small q. The appearance of the oscillations can 
be traced back to the calculation of the pseudized dis- 
paced valence electronic density which has a discontin- 
uous second derivative at a matching radius. However, 
these oscillations do not influence the final OF-AIMD re- 
sults because they appear for g-values bigger than those 
corresponding to the -Ecut- 
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FIG. 1: Non-coulombic part of the pseudopotential for Al at 
T = 943 K. The continuous hue is the OFDFT-PS used in 
the OF-AIMD simulations, while the dashed line is stands for 
the LRT-PS used for the CMD simulations. 



B. Static properties. 

The static structure factors, S(q), obtained from the 
simulations are shown in Figure 0, which also shows the. 
corresponding experi mant al data measured by neutronEHI 
and X-ray difFractiontilEa experiments. The experimen- 
tal data show small differences in the region 2 
< g < 5 A~^, with the neutron values being slightly big- 
ger that the X-ray ones, whereas the OF-AIMD results 
stand remarkably well between both sets, although some- 
what closer to the X-ray data. The insets of the figures 
show that the OF-AIMD results in the small g-region are 
also in good agreement with the experimental X-ray re- 
sults. The figures also include the S{q), obtained from 
the CMD simulations performed with the interatomic 
pair potential derivjed from the LRT-PS and from the 
QHNC method.ESy Although these calculated S{q) rea- 
sonably reproduce the experimental data the agreement 
with experiment is much better for OF-AIMD. 

Extrapolation of S{q) to q ^ allows the isother- 
mal compresibility, kt, to be estimated from the relation 
5(0) — phsTnT- A least squares fit of S{q) — sq + S2q^ 
to the calculated S{q) for q- values up to 1 yields the 
result kt,of-aimd = 2.37 (in 10~^^ m^ Nw~^ unitsi for 
T = 943 K, which is close to the experimental valueE3 kt 
= 2.43. In contrast, both the LRT-PS and QHNC inte- 
rionic pair potentials lead to much higher values, namely 
kt,lrt-cmd = 6.5 and kt,qhnc-Cmd = 7.4 respectively. 

The ionic and electronic static structure of liquid 
Al near melting has also been calculated by Anta et 




q(A-^) 



FIG. 2: Static structure factors of liquid Al at (a) 943 K 
and (b) li3^-i K. Full circles: experimental X-ray diffrac- 
tion datatEjI'LJ Open circles; experimental neutron diffrac- 
tion data.E3 Continuous line: OF-AIMD simulations. Dashed 
lines: LRT-CMD simulations. Dotted lines: CMD simula- 
tions with the QHNC potential. The insets show the low-g 
behaviour. 

using the OF-AIMD method with a kinetic en- 
ergy functional which describes the corxect linear and 
quadratic response of the electron gasCS and a local 
ionic pseudopotqnJtiaJ constructed from a non-local ionic 
pseudopotential. E3'c£l Their results for the static structure 
factor closely followed the experimental one. 

C. Dynamic properties. 

1. Single-particle dynamics. 

The most complete information about the single- 
particle properties is provided by the self-intermediate 
scattering function, Fs{q,t), which probes the single- 
particle dynamics over different length scales, ranging 
from the hydrodynamic limit {q —>■ 0) to the free-particle 
limit {q oo). In the present simulations, this magni- 
tude has been obtained by 



F,{q,t) = — e-'QRAt+ta)^^qRAto)^ (26) 
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t(ps) 

FIG. 3: Self intermediate scattering functions, Fs{q,t), at 
several g-values, for liquid aluminium. Continuous, dashed 
and dotted lines: OF-AIMD, LRT-CMD and QHNC-CMD 
simulations respectively at T= 943 K. Dash-dotted line: OF- 
AIMD simulations at T=1323 K 

and in figure || we show the results obtained for sev- 
eral g-values at T= 943 K and 1343 K. It shows the 
typical monotonic decrease with time; moreover, the re- 
sults are very similar to those of the LRT-CMD and 
QHNC-CMD simulations, although the latter ones show 
a slightly slower decay with time. An increase in tem- 
perature leads to increased rate of decay. 

Closely related to the Fs{q,t) is the velocity autocor- 
relation function (VACF) of a tagged ion in the fluid, 
Z(t), which can be obtained as the q — > limit of the 
first-order memory function of the Fs{q,t). However, in 
the present simulations it is more easily obtained from 
its definition 

Zit)^{Mt)MO))/{vl) (27) 

which stands for the normalized VACF. The results are 
shown in Figure ^ along with those derived fron the LRT- 
CMD and QHNC-CMD simulations. The results dis- 
play the typical backscattering behaviour, which is more 
marked for the OF-AIMD simulations, but the results 
of the three simulations are rather similar. The main 
features of the obtained Z{t) are comparable |-to-4iiosa 
obtained for other simple metals near mclting,E3c3e3'c3 
namely: (i) a first minimum about 0.20 deep and (ii) a 
rather weak following maximum peaking close to cero. 




t(ps) 



FIG. 4: Normalized velocity autocorrelation functions. Con- 
tinuous, dashed and dotted lines: OF-AIMD, LRT-CMD and 
QHNC-CMD simulations respectively at T= 943 K. Dash- 
dotted line: OF-AIMD results for T= 1323 K 

The self-diffusion coefficient, D, is readily obtained from 
either the time integral of Z{t) or from the slope of the 
mean square displacement SR'^{t) = (|^i(<) — ^i(O)p) of 
a tagged ion in the fluid, as follows 

1 f°° 

/ Z{t)dt; D ^ Mm 6R^{t)/m (28) 

/3m Jo t-oo 

and the results for D are given in Table ||. The two 
routes for D lead to practically the same value, namely 
-DoF-AiMD — 0.49 A^/ps; which is somewhat smaller 
than the mean value of 0.55 A^/ps obtained iii a pre- 
vious OF-AIMD calculation with 205 particles.El Unfor- 
tunately, to our knowledge, no experimental results are 
yet available for the diffusion coefficients of liquid Al 
at any thermodynamic state. However, we cap-, com- 
pare with the results of a KS-DFT calculationEl per- 
formed for liquid Al near the tripe point, using 64 par- 
ticles and a nonJpcai Bachelet-Hamann-Schluter type 
pseudopotentiahElSEj this calculation lead to a value 
-Dks-dft — 0.60 A^/ps derived from the slope of the 
corresponding mean square-displacement. Recently, an- 
other KS-DFT calculatioiJd for liquid Al at 1000 K, 
using 64 particles and ultrasoft Vanderbilt pseudopo- 
tentials gave for D values within the range 0.52-0.68 
A^/ps, derived from the slope of the mean square dis- 
placement. Our OF-AIMD simulations, with a small 
number of particles and/or a small number of configura- 
tions, suggest that the self-diffusion coefficients obtained 
from the 5R^{t) tend to be greater than those obtained 
by integration of the Z(i), and as the number of parti- 
cles and/or configurations is increased, the value for the 
self-diffusion coefficient is decreased. More extensive KS- 
DFT simulations would probably lead to a smaller value 
of D closer to that obtained in the present OF-AIMD 
simulations. The values obtained from LRT-CMD (0.58 
AVps ) and QHNC-CMD (0.55 AVps ) simulations are 
also rather similar and slightly bigger than the OF-AIMD 
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FIG. 5: (a) Normalized half- width of Ss{q,i^), relative to its 
value at the hydrodynamic limit, for liquid aluminium at T 
= 943 K (continuous line) and 1323 K (dash-dotted line). 
The dashed line stands for the free-particle limit, (b) Same 
as before, but for the normalized peak value Ss{q,(^ = 0), 
relative to its value at the hydrodynamic limit 

result. The CMD simulations of Ebbsjo et azB using sev- 
eral pair potentials gave values within the range 0.41-0.45 

A^/ps. 

By Fourier transforming the Fs{q,t) we obtain the 
self-dynamic structure factor, Ss{q,tij), which for all q- 
values, exhibits a monotonic decay with frequency, from 
a peak value at oj = 0. Ss{q-,Lo) can be characterized 
by the peak value S's(g,ti; = 0), and the half-width at 
the half- maximum, ^1/2(9)- These parameters are usu- 
ally reported normalized with respect to the values of 
the hydrodynamic (g ^ ) limit, by introducing the 
dimensionless quantities S((z) = Trq^DSs{q,uj = 0) and 
A(g) = uJi/2(<l)l<f'D, where Wi/2('?)/q^ can be inter- 
preted as an effective g-dependent diffusion coefficient 
D{q). For a liquid near the triple point, A((7) usually 
exhibits an oscillatory behaviour whereas, in a dense 
gas it decreases monotonically from unity at g = to 
the 1/q behaviour at large q. Fig. || shows the OF- 
AIMD results for /S.{q) and the corresponding re- 

sults from the LRT-CMD and QHNC-CMD simulations 
are rather similar and are not shown. The results for 
Aof-aimd(<z), show that for both temperatures, the hy- 
drodynamic limit is reached from below, with a mini- 
mum at around g ~ Qp, followed by a maximum and 
by a gradual transition, for greater g-values, to the free- 
particle limit. Notice that for the higher temperature, 
the oscillations are heavily damped and the free particle 



limit is approached quickly. This oscillating behaviour 
of Aof-aimd('?) for small and intermediate q- values has 
been reported by several authors and has been attributed 
to the coupling of th e ,, sipgloj part icle motion to other 
modes in the system.c3E3E3£Zl'E3'E2l On the other hand, 
the results for S(fc) reflect greater sensitivity to changes 
in temperature, with the diffusive limit reached from be- 
low for T=943 K and from above for T=1323 K. We note 
that similar features to those obtained in this paper were 
obtained earlier by Torcini et aZ.Ej in their CMD study 
of liquid lithium near melting usipa the interatomic pair 
potential proposed by Price et aZ.Ell 



2. Collective dynamics. 

The intermediate scattering function, F{q,t), embod- 
ies the information concerning the collective dynamics of 
density fluctuations over both the length and time scales. 
It is defined as 



F{q,t) 




^-iqRjit+to) 




(29) 



and in figures we show, the results from the present 
OF-AIMD simulations for several g- values. F{q, t) ex- 
hibits oscillatory behaviour which persists up to g w 
3(7j,/5, with the amplitude of the oscillations being 
stronger for the smaller g-values. This is typical be- 
haviour found for other simple liqui d j p a ta ls near melt-, 
ing, by either computer simulationg3^c3 or theory.o 
Different behaviour is seen for the results in the same 
(j-range obtained from the LRT-CMD and QHNC-CMD 
simulations, with F(q, t)'s whose contact values, given by 
F{q,t = 0) = S{q), are more than double and, more im- 
portant, with a diffusive component playing a domiaant 
role. The corresponding MD results of Ebbjso et alB for 
the F{q,tys have better contact values but also display 
an important diffusive component. 

Closely connected to the F{q, t) is the dynamic struc- 
ture factor, S(q,uj), which is obtained by a time Fourier 
transform of the F{q, t) (with an appropriate window to 
smooth out truncation effects). Its importance lies in 
its direct connection to the inelastic neutron scattering 
or the IXS data. The results obtained for the S{q^ lu) 
are shown in figures for a range of wavevectors up 
to w 2.5qp. The dynamic structure factor shows well 
defined sidepeaks, indicative of collective density excita- 
tions, up to q « 1.6 which amounts to « 3qp/5. The 
results qualitatively reproduce the shape of the experi- 
mental IXS dataO with some small discrepancies in the 
heigths and positions of the peaks. Similar results, but 
with a better description of the central peak at the low- 
est g-values, were also obtained in the CMD simulations 
of Ebbsjo et alQ However, it must be stressed that their 
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FIG. 6: Normalized intermediate scattering functions, F{q,t), 
at several g-values, for liquid aluminium at T = 943 K. Con- 
tinuous line: OF-AIMD simulations. Dashed line: LRT-CMD 
results. Dotted line: QHNC-CMD results. 



F{q,tys were previously fitted to an analytical expres- 
sion interpolating among the ideal gas, viscoelastic and 
hydrodynamic models, and therefrom the corresponding 
S{q,uj) were derived. The strongly diffusive character 
of the F{q,tys obtained from both the LRT-CMD and 
QHNC-CMD simulations, give rise to S{q,uj) which de- 
cay rather quickly, with hardly discernable side peaks. 
This is because the side peaks are located at smaller po- 
sitions, given by qcs{q), where Cs{q) is the generalized 
adiabatic sound velocity (see below), which is too small 
because of the large values of S{q) at those q- values. 

From the positions of the sidepeaks, LJm{q), the dis- 
persion relation of the density fluctuations has been ob- 
tained and this is shown in figure |l^ for the state at 
T = 943 K, along with u)i{q), which is the dispersion 
relation obtained from the maxima of the longitudinal 
current correlation function, Ji(q,uj) = u>'^S{q,u!). Note 
that in the hydrodynamic region (small q), the slope of 
the dispersion relation curve is the adiabatic sound ve- 
locity, Cs{q) = vth\/l/S{q), with vth = {l3m)-^/'^ being 
the thermal velocity and 7 is the ratio of the specific 
heats. In the limit q ^ 0, Cs{q) reduces to the bulk 
adiabatic sound velocity and determines the slope of the 
dispersion at g — > 0. By extrapolating the OF-jAIMD 
results for S{q) and using the experimental valueE3 of 7 
w 1.25, we obtain a value of « 4850 m/s for the bulk 
adiabatic sound velocity whieh compares reasonably well 
with the experimental valueca of w 4700 m/s, near the 
triple point. Figure [IQ shows a positive dispersion, i.e. 



FIG. 7: Same as the previous figure 



an increase of toi (q) with respect to the values predicted 
by the hydrodynamic adiabatic speed of sound, with a 
maximum located around 0.4 A~^. Similar behaviour 
has also been obtained by Scopigno et a/.Lj from their 
experimental IXS results for liquid Al at T=1000 K and 
ha S| b|jf^ : ^ bserved in other liquid metals: Rb, Cs, Li and 

Another interesting dynamical magnitude is the trans- 
verse current time correlation function, Jt{q,t), which 
is not associated with any measurable quantity and can 
only be determined by means of MD simulations. It pro- 
vides information on the shear modes and is defined as 



Jtiq,t) = - {j;iq,0)jAq,t)) 



(30) 



where jx{q,t) — Sj^^i^JlO^ ^^^^jW is the transverse 
current. The shape of Jt{q, t) evolves from a gaussian, in 
both q and t, for the free particle g — > 00 limit, towards a 
gaussian in q and exponential in t for the hydrodynamic 
limit (q — > 0), i.e. 



Jt{q 



Q t) = J-e-''"!*!/"" 



(31) 



where 77 is the shear viscosity. For intermediate g-values, 
Jt{q,t) exhibits a more complicated behaviour, as shown 
in Figure |ll| where OF-AIMD results for liquid Al near 
melting are shown. Note that for the smallest g- value 
reached by the simulation: q — 0.29 A~^, the corre- 
sponding Jt{q, t) takes on negative values, which by eqn. 
( ^ ) means that it is already beyond the hydrodynamic 
regime. The associatted spectrum, Jt{q,uj), plotted in 
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FIG. 8: Dynamic structure factor, S{q,ui), for several q- 
values, for liquid aluminium at T = 943 K. Continuous line: 
OF-AIMD simulations. Dashed line: LRT-CMD results.r*ot- 
ted line: QHNC-CMD results. Full circles: Experiment^. 



Figure |l^, shows an inelastic peak which already exists 
at q = 0.29 « O.llgp ; as q increases the peak be- 
comes better defined and it persists to q values around 
3qp, although it has already disappeared for the largest 
g-value considered. Note that the associatted peak fre- 
quency increases with g up a maximum value at g ~ g^, 
and then flattens at larger q as Jt{q,uj) evolves towards 
a gaussian shape. This behaviour closely parallels that 
observed for the alkali-metals where the inelastic peak 
appears for q > 0.07(7p.c3 

Similar results are also obtained by the LRT-CMD ap- 
proach, but Jt{q,t) decays more slowly and the minima 
are less marked. This leads to a spectrum Jt{q,uj) where 
the peaks are less marked and, in fact, there is no peak 
at g = 0.29 A^i. 

From the results for Jt{q,t) we can rpariihL obtain the 
shear viscosity coefficient, rj as follows Jl3l£3'E£l The mem- 
ory function representation of the Jt{q,t): 



FIG. 9: Same as the previous figure 
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FIG. 10: Dispersion relation for liquid Al at T =943 K. Open 
triangles: peak positions, tJm(g), from the OF-AIMD S{q,u}). 
Open circles: peak positions, LJi{q), from the maxima of the 
OF-AIMD longitudinal current, J; (g, oj). Full line: Linear dis- 
persion with the hydrodynamic sound velocity, v=4700 m/s. 



pm 



z H fi{q,z) 

pm 



(32) 



where the tilde denotes the Laplace transform, introduces 
a generalized shear viscosity coefficient, f]{q, z). The area 
under the normalized Jt{q, t), gives (3m Jt{q, z — 0) from 
which values for z = 0) can be obtained which when 
extrapolated to q — give the usual shear viscosity co- 
efficient, rj. Results for ij presented in Table || compare 
favourably with the available experimental data.Ej For 
comparison—we note that the KS-DFT simulations of Alfe 
and GillanEJ gave values in the range 1.4 - 2.2 GPa • ps. 



IV. CONCLUSIONS. 

Several dynamic properties of liquid aluminium have 
been calculated at two thermodynamic states close to the 
triple point. The simulations have been performed us- 
ing the orbital free ab-initio molecular dynamics method, 
showing the feasibility of this technique to calculate sev- 
eral time correlation functions, allowing a comprehensive 
study of the dynamical properties. Furthermore, agree- 
ment with the available experimental data is quite satis- 
factory. 
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FIG. 11: OF-AIMF transverse current correlation function, 
Jt{q,t), at several g-values (in A~^) for liquid Al at T =943 
K. (a) g=0.29 (full curve), g=0.42 (dashed curve), g=0.98 
(dotted curve), q= 1.45 (dash-dotted curve), (b) 5=2.70 (full 
curve), q= 3.80 (dashed curve), q= 5.54 (dotted curve), q= 
9.02 (dot-dashed curve). 
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FIG. 12: OF-AIMF transverse current correlation spectrum, 
Jt{q,oj), at several g- values (in A~^) for liquid Al at T =943 
K. (a) g=0.29 (fuU curve), q=QA2 (dashed curve), g=0.98 
(dotted curve), q= 1.99 (dash-dotted curve), (b) 5=2.70 (full 
curve), g= 3.80 (dashed curve), q= 5.54 (dotted curve), q= 
7.03 (dot-dashed curve), q= 9.02 (pluses). 



We have also presented a method for producing, from 
first principles, local pseudopotentials for use with the 
orbital free density functionals. While the ultimate goal 
of the method would be to use the atomic number of 
the atoms as the only input data, this has not yet been 
achieved as the present calculations also require the ex- 
perimental number density of the system for calculating 
the local pseudopotential and for performing the simula- 
tions. However, we stress that starting from very basic 
information, the present scheme allows the determination 
of the static and dynamic properties of the system. 

Finally, we emphasize that in the present scheme, the 
calculation of the pseudopotential is coupled to the par- 
ticular functional adopted for the total potential energy 
of the system. This means that different kinetic en- 
ergy functionals would lead to different pseudopotentials. 
Consequently, this field is open to further improvements 
in the description of the kinetic energy and, therefore, 
also in the corresponding local pseudopotential. 



APPENDIX A: THE KINETIC ENERGY 
FUNCTIONAL 

We consider the kinetic energy functional 



TM^Tw[p]+Tp[p] 



where 



3 

To 



dr p{r 



,5/3-2/3 



~k{r) 



(Al) 



(A2) 



k{r) = {2k' 



\3 



dsk{s) ujfj{2k°p\r- s\) = k{r) * uj0{2k°pr) 

(A3) 



k{r) = (37r2)i/3p(j;-)/3 



(A4) 



In the limit of small deviations from a uniform sys- 
tem, we wish to recover the LRT result. Equating the 
Fourier transform (FT) of the second functional deriva- 
tive of Ts[p] with respect to p{r) for p{r) = po, to the 
inverse of the Lindhard response function, gives for the 
weight function 



6/3' - f /3 + y) + 4/3 ("^ - 2/3 ) U^{rj) + 20'ZJp{^Y 



^(lA,(ry)-3ry^) 
where rj = q/2kp, Tjp is the FT of ojp and 



(A5) 



7rL(??) 



2V 



In 



1-77 



(A6) 



is the noninteracting hom oge neous electron gas response 
function. Taking in eqn. ( [Aq ) the solution which satisfies 
the normalization condition ujfj{ri = 0) = 1, and with /3 
within the range < /3 < 5/6 so that the power of p{r) 
in eqn. (A2) is positive, the weight function is given by 
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and in terms of the modified weight function 



5 

3^ 



^ ^'(S -3/3)2 + 5(71^1 (77) -1 



3/3 



3772). 
(A7) 



Requiring tZJ^ to be real places a stricter limit on (3: 
f3 < 0.5991. With this choice of weight function, the 
functional recovers the LRT limit, and in the limit of uni- 
form density it reduces to the Thomas-Fermi functional. 
In the limit 77 —> 00 we have 



(A8) 



where 



C. = 2-A + _Lyi73^5^ (A9) 

The constant Ci gives rise to a Dirac delta function 
in the real space; therefore it is convenient to define a 
"modified" weight function 



'^f^i'n) = LUfsifj) - Ci 



(AlO) 



so that every convolution involving ujp, such as in eqn. 
([A^), becomes 



G(f) * ujf3{2k%r) = CiG(f) -I- G(f) * w^(24r) . (All) 

An important limit is when the mean electron den- 
sity, and therefore kPp , vanishes as for instance in a finite 
system. Now, the convolutions involving the "modified" 
weight function vanish because 77 = q/2kp 00 and uj^rj) 
vanishes. Consequently, k{r) = Ci k{f), and the kinetic 
energy functional becomes Ts[p] — Tw[p] + Ci TTF[p]^ 
and when /5 = 4/9, Ci = 0. 



APPENDIX B: THE POSITION DEPENDENT 
CHEMICAL POTENTIAL. 

The functional derivative of eqn. (^ gives 



H{f) = fj,y^i/) + npi/) V^^t{r) + VH{r) + V^c[p{r)] 

(Bl) 

where 



p.w{r) 



Vp{r) 



1 V'^p{t^ 
8 p{rf 4 p{r) 



Vh{t^ = J ds 



\r-s\ ' 



(B2) 



(B3) 



M/3(^) 
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(5/3-2/3) p(r)2/3-2/3^(^)2 ^ 
2/3(37r2)i/3p(f)/3-i/^(f) ] 



with 



h{r) = f{r)*G:p{2k'j,r) 



(B4) 



(B5) 



where f{r) — k{r) p(r)^/3-2/3 ^ rpj^^ product ^{r)'ip{r), is 
the "driving force" for the dynamical minimization of the 
energy functional, see eqn. dl3) . If the various powers of 
the density appearing in ^,p{ryip{f) are to remain positive 
so that this driving force is not to diverge in regions where 
the density vanishes, then 1/2 </3< 7/12. In practice 
we have found that for /3 = 0.51 the minimization has 
always proved possible. 



APPENDIX C: CONSTRUCTING THE LOCAL 
PSEUDOPOTENTIAL FROM AN INFINITE 
SYSTEM 

In an infinite system most of the previous expressions 
diverge, because the integrals extend to all space, and the 
integrands don't vanish for large distances. Moreover, 
the normalization constraint must be redefined, because 
the total number of electrons is infinite. 

To avoid these problems, one has to take into account 
the similar divergencies which appear in the "ionic" part 
of the total energy. This amounts to use the "difference" 
functions, which are obtained by substracting from the 
total functions their corresponding limits for large dis- 
tances, and redefining the chemical potential as the La- 
grange multiplier associated to the normalization of the 
"displaced" density n{r) — p{r) — po- In this way we 
define the following functions: 

. A„(f) = p(f)" -p^ = (po + n{r)r - Po 
In particular, for a = 1, Ai(r') = n{r), the dis- 
placed density. 

. x(r) = k{r) - (3^2)1/3^/3 

. </)(f) = fir) (37r2)i/3pV3-^* 

.r72(r') = /7(f)-(l-Ci)(3^2)i/3^5/3-/5 

.~^,{r)^^,,{r)-\{k^pf 

• Vc^t{r) = Vc^t{r) - vic\\{r), where v-;^c\\{r) is the po- 
tential created by a uniform background of positive 
charge with density po- 



Vxcir) = Vxc[po + 77 (r*)] - Vxc[po] 
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In terms of these functions, the Euler equation now In the last equation, 772(7^ = ipir) *Cjp{2k^pr) and 
becomes: 



Vc^t{r) + / dsn{s)/ \r-s\+ v^dr) + 



Vn(f) 



(po + nif)y 



1 V^n(r) 
4 po + n{r) 



+ fipir) - = (CI) 



with 



A/3 (^) = T7T [ Aa (r^O + MS (r ) + Ac (r ) + (r ) ] , (C2) 
where 
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AA(r) = U - 2/3 ) p^-'' \2i3nY^^p^„m + xir) 



(C3) 



AB(rO = - 2/3") A2/3_2/3(r) [(37r2)^/V(? + xir) ' 



(C4) 



</.(r) = pT-'^'m + {3ny/'p'o^,/3-2p{f^ 

+ A5/3~2M0X(^ (C7) 



Summarizing, to evaluate p^[3{r) from the displaced 
density n(f) we take the following steps: 

1. Compute Ac,(f) for a = 5/3-2/3,2/3-2^,2/3-/3, 
and /3 — 1 

2. Compute x(r) and FT to obtain xi?) 

3. Compute x(g) = Cix(g) + x{q)^i3{ql'2k%) and by 
inverse FT obtain xl?^) 

4. Compute (/)(r) according to eqn. ( |C7| ) and FT to 
obtain ^((t) 

5. Compute ri2{q) ~ (f){q)ujp{q/2kp) and inverse FT 
to obtain 772 (r) 



6. Apply eqns. ( |C2| - |C6[) to obtain ppir) 



When our system is 1 atom in a jellium-vacancy the 

- /-^ fo 2\i/3r 2/3-/3-/^N , fr, 2\i/3 P\ external potential is given by 

/ic (r) = 2/3Ci(37r ) /•'[po xW + (Stt ) '•'p^A2/3-/3W ^ b J 



and 



+ A2/3-/3(r)x(f)] , (C5) 



poir) - 2/3(1 - Ci)(37r2)2/3p5/3-2/3^^_^(.^ 



+2/3(37r2)i/3[po + n(f/-S2(rO ■ (C6) 



i.e., the sum of the potential created by the cavity and 
the ionic pseudopotential. Substituting into eqn. (CI) 
we directly obtain Vpsif) once we compute all the other 
terms, which are calculated from n{r). 

Note that all the functions have spherical symmetry, 
which leads to simple expressions for the gradient, the 
laplacian and also the Fourier tranforms. 
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TABLE I: Thermodynamic states studied in this work, along 
with some simulation details. 





T(K) 


N 


p (A-«) 


£cut(Ryd) 


OF-AIMD 


943 


500 


0.05290 


30.25 


OF-AIMD 


1323 


500 


0.05071 


29.25 


LRT-CMD 


943 


600 


0.05290 




LRT-CMD 


1323 


500 


0.05071 




QHNC-CMD 


933 


800 


0.05331 





TABLE II: Isothermal compresibility kt (in 10"" m^ N"^), 
self-diffusion coefficient D (in A^/ps) and shear viscosity coef- 
ficient r] (in GPa-ps) of liquid Al at the thermodynamic states 
studied in this work. 





T(K) 




D 


V 


OF-AIMD 


943 


2.37 


0.49 


1.38 


OF-AIMD 


1323 


2.38 


1.05 


0.85 


LRT-CMD 


943 


6.57 


0.58 


1.24 


LRT-CMD 


1323 


6.32 


1.14 




QHNC-CMD 


933 


7.45 


0.55 


1.36 


Experiment 


933 


2.43 " 




1.26 
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